exposure_name <- c("met_homovanillic_acid", "met_histidine", "met_hydroxyphenyllactic_acid")
pfas_c <- c("pfas_pfda", "pfas_pfhpa", "pfas_pfna", "pfas_pfuna", "pfas_pfoa")
pfas_s <- c("pfas_pfhps", "pfas_pfhxs", "pfas_pfos")
outcome <- c("mic", "daystomic", "daystomac", "rapid", "daystorapid","daystohyp")
table1::table1(~factor(mic) + factor(rapid), full_data %>% select(exposure_name, pfas_c, pfas_s, outcome))
| Overall (N=40) |
|
|---|---|
| factor(mic) | |
| 0 | 30 (75.0%) |
| 1 | 10 (25.0%) |
| factor(rapid) | |
| 0 | 29 (72.5%) |
| 1 | 11 (27.5%) |
table1::table1(~., full_data %>% select(exposure_name, pfas_c, pfas_s, outcome[c(2,3,5)]))
| Overall (N=40) |
|
|---|---|
| met_homovanillic_acid | |
| Mean (SD) | 0.100 (0) |
| Median [Min, Max] | 0.100 [0.100, 0.100] |
| met_histidine | |
| Mean (SD) | 66.8 (7.42) |
| Median [Min, Max] | 65.9 [51.4, 82.0] |
| met_hydroxyphenyllactic_acid | |
| Mean (SD) | 0.559 (0.503) |
| Median [Min, Max] | 0.374 [0.100, 1.68] |
| pfas_pfda | |
| Mean (SD) | 0.252 (0.107) |
| Median [Min, Max] | 0.240 [0.141, 0.530] |
| pfas_pfhpa | |
| Mean (SD) | 0.299 (0.247) |
| Median [Min, Max] | 0.225 [0.0707, 1.10] |
| pfas_pfna | |
| Mean (SD) | 1.17 (0.645) |
| Median [Min, Max] | 0.965 [0.290, 2.80] |
| pfas_pfuna | |
| Mean (SD) | 0.155 (0.128) |
| Median [Min, Max] | 0.120 [0.0707, 0.790] |
| pfas_pfoa | |
| Mean (SD) | 3.16 (1.38) |
| Median [Min, Max] | 2.90 [0.960, 9.30] |
| pfas_pfhps | |
| Mean (SD) | 0.273 (0.124) |
| Median [Min, Max] | 0.245 [0.141, 0.710] |
| pfas_pfhxs | |
| Mean (SD) | 3.37 (3.04) |
| Median [Min, Max] | 2.45 [0.890, 15.0] |
| pfas_pfos | |
| Mean (SD) | 9.13 (3.79) |
| Median [Min, Max] | 8.20 [3.10, 22.0] |
| daystomic | |
| Mean (SD) | 3110 (1670) |
| Median [Min, Max] | 3650 [298, 5510] |
| daystomac | |
| Mean (SD) | 3470 (1500) |
| Median [Min, Max] | 4000 [629, 5510] |
| daystorapid | |
| Mean (SD) | 3130 (1560) |
| Median [Min, Max] | 3290 [348, 5510] |
G = full_data %>% select(exposure_name[-1]) %>% as.matrix() %>% scale()
Z = list(Carboxylic = full_data %>%
select(sample_id, all_of(pfas_c)) %>%
column_to_rownames("sample_id") %>%
as.matrix() |> scale() ,
Sulfonic = full_data %>%
select(sample_id, all_of(pfas_s)) %>%
column_to_rownames("sample_id") %>%
as.matrix() |> scale() )
Y = full_data %>%
dplyr::select(all_of(outcome)) %>%
mutate_at(.vars = c(outcome),
.funs = ~unlist(.) %>% as.numeric()) %>%
mutate_at(.vars = c(outcome[c(2,3,5,6)]),
.funs = ~scale(.))
df <- data.frame()
plot <- vector("list", length(length(outcome)))
for (i in c(2,3,5,6)){
fit <- LUCIDusM::est_lucid(
G = G,
Z = Z,
Y = Y[,i],
K = rep(2, length(Z)),
family = "gaussian",
max_itr = 200,
useY = TRUE,
modelNames = rep("EEV", length(Z)))
fit_reordered <- reorder_lucidM(fit, reference = c(1,2))
# fit_reordered <- fit
est <- data.frame(
outcome = colnames(Y)[i],
Exposure1_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,2],
Exposure1_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,2],
Exposure2_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,3],
Exposure2_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,3],
# Exposure3_to_Carboxylic =
# fit_reordered$res_Beta$Beta[[1]][,4],
# Exposure3_to_Sulfonic =
# fit_reordered$res_Beta$Beta[[2]][,4],
Carboxylic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[2],
Sulfonic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[3])
df <- df %>% bind_rows(est)
(plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
sankey_colors =
sankey_colors,
text_size = 12))
names(plot)[i] <- paste0(colnames(Y)[i])
htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}
for (i in c(1)){
fit <- LUCIDusM::est_lucid(
G = G,
Z = Z,
Y = Y[,i],
K = rep(2, length(Z)),
family = "binomial",
max_itr = 200,
useY = TRUE,
modelNames = rep("EEV", length(Z)))
fit_reordered <- reorder_lucidM(fit, reference = c(1,1))
# fit_reordered <- fit
est <- data.frame(
outcome = colnames(Y)[i],
Exposure1_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,2],
Exposure1_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,2],
Exposure2_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,3],
Exposure2_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,3],
# Exposure3_to_Carboxylic =
# fit_reordered$res_Beta$Beta[[1]][,4],
# Exposure3_to_Sulfonic =
# fit_reordered$res_Beta$Beta[[2]][,4],
Carboxylic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[2],
Sulfonic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[3])
df <- df %>% bind_rows(est)
(plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
sankey_colors =
sankey_colors,
text_size = 12))
names(plot)[i] <- paste0(colnames(Y)[i])
htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}
for (i in c(4)){
fit <- LUCIDusM::est_lucid(
G = G,
Z = Z,
Y = Y[,i],
K = rep(2, length(Z)),
family = "binomial",
max_itr = 200,
useY = TRUE,
modelNames = rep("EEV", length(Z)))
fit_reordered <- reorder_lucidM(fit, reference = c(2,2))
# fit_reordered <- fit
est <- data.frame(
outcome = colnames(Y)[i],
Exposure1_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,2],
Exposure1_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,2],
Exposure2_to_Carboxylic =
fit_reordered$res_Beta$Beta[[1]][,3],
Exposure2_to_Sulfonic =
fit_reordered$res_Beta$Beta[[2]][,3],
# Exposure3_to_Carboxylic =
# fit_reordered$res_Beta$Beta[[1]][,4],
# Exposure3_to_Sulfonic =
# fit_reordered$res_Beta$Beta[[2]][,4],
Carboxylic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[2],
Sulfonic_to_outcome =
fit_reordered$res_Delta$fit$coefficient[3])
df <- df %>% bind_rows(est)
(plot[[i]] <- plot_lucid_in_parallel_plotly(fit_reordered,
sankey_colors =
sankey_colors,
text_size = 12))
names(plot)[i] <- paste0(colnames(Y)[i])
htmlwidgets::saveWidget(plot[[i]], file = fs::path(dir_figure,str_c("lucid_in_parallel_", colnames(Y)[i],".html")))
}
writexl::write_xlsx(df, fs::path(dir_results, "LUCID/lucid_in_parellel_est.xlsx"))
plot[[1]]
plot[[2]]
plot[[3]]
plot[[4]]
plot[[5]]
plot[[6]]
G = full_data %>% select(exposure_name[-1], "met_glycine", "met_2_methyl_3_hydroxybutyric_acid") %>% as.matrix() %>% scale()
fit <- LUCIDusM::est_lucid(
G = G,
Z = Z,
Y = Y[,4],
K = rep(2, length(Z)),
family = "gaussian",
max_itr = 200,
useY = TRUE,
modelNames = rep("EEV", length(Z)))
fit_reordered <- reorder_lucidM(fit, reference = c(2,2))
p <- plot_lucid_in_parallel_plotly(fit_reordered,
sankey_colors =
sankey_colors,
text_size = 12)
htmlwidgets::saveWidget(p, file = fs::path(dir_figure,"lucid_in_parallel_rapid_with_4met.html"))
p